# aem.two_d.transient.base_aquifer
#
# Copyright (c) 2000, V.Kelson and M.Bakker
#           (c) 2001, Willem Jan Zaadnoordijk

import geometry
SearchableList = geometry.two_d.collections.SearchableList

import aem.container
RowList = aem.container.RowList
ColumnList = aem.container.ColumnList

import base
from aem import matrix

class Aquifer(base.Element):
	"""Aquifer base class for single aquifer

	This is the base class for transient AquiferSystem objects.
	It builds the aquifer system in terms of potentials only;
	no aquifer properties are provided.
	This is not quite true: the aquifer diffusivity is needed
	(it provides the necessary relation between the spatial
	and the temporal dimensions and) it has a default value 1
	Although instances of this class may be created,
	the typical application will override the constructor to
	add aquifer properties, uniform flow, and other features.
	"""
	def __init__(self,kH_S=1.0,C=0.0,parent=None):
		"""Constructor

		Accepts the constant (C) for the potential.
		"""

		self.kH_S=kH_S        # aquifer diffusivity kH/S
		self.initialSteadyState=aem.two_d.single.base_aquifer.Aquifer(C=C)
		self.elementList = SearchableList([self])
		base.Element.__init__(self,parent)
		self.givenElement = 0
	def addElement(self,el):
		self.elementList.append(el)
		# add check to see if element el is based on transient element
		# if not then in stead of above statement:
		# self.initialSteadyState.elementList.append(el)
	def potential(self,pt,t=0):
		rv = 0.0
		for e in self.initialSteadyState.elementList: rv = rv + e.potentialContribution(pt)
		for e in self.elementList: rv = rv + e.potentialContribution(pt,t)
		return rv
	def discharge(self,pt,t=0):
		rv = geometry.XY(0.0,0.0)
		for e in self.initialSteadyState.elementList: rv = rv + e.dischargeContribution(pt)
		for e in self.elementList: rv = rv + e.dischargeContribution(pt,t)
		return rv
	def recharge(self,pt,t=0):
		rv = 0.0
		for e in self.initialSteadyState.elementList: rv = rv + e.rechargeContribution(pt)
		for e in self.elementList: rv = rv + e.rechargeContribution(pt,t)
		return rv
	def extraction(self,t=0):
		rv = 0.0
		for e in self.initialSteadyState.elementList: rv = rv + e.extractionContribution()
		for e in self.elementList: rv = rv + e.extractionContribution(t)
		return rv
	def head(self,pt,t=0):
		return self.initialSteadyState.potentialToHead(pt,self.potential(pt,t))
	def aquiferDiffusivity(self,pt,t):
		return self.kH_S

	def solveIteratively(self,factor=1.0,recursive=0,t=0):
		'''Compute solution iteratively'''
		self.solve(factor,[self])
		for e in self.elementList[1:len(self.elementList)]:
			if recursive>1 and len(e.elementList)>1:
				e.solveIteratively(factor,recursive-1,t)
			else:
				e.solve(factor,e.elementList,t)

	# Overload these in derived classes
	def headToPotential(self,pt,head):
		"""Returns head; must be overloaded"""
		return head
	def potentialToHead(self,pt,pot):
		"""Returns potential; must be overloaded"""
		return pot
	def velocity(self,xyzpt,t=0):
		'''Returns velocity; must be overloaded'''
		velo = geometry.XYZ(0.0,0.0,0.0)
		return velo

